home *** CD-ROM | disk | FTP | other *** search
/ Mac Mania 4 / MacMania 4.toast / / Demo's / Igor Demo Pro / 1 PutContentsIn Igor Pro Folder / WaveMetrics Procedures / Math Functions / AreaXY next >
Text File  |  1996-04-01  |  2KB  |  66 lines
  1. #include <BinarySearch>
  2.  
  3. | AreaXY 1.2
  4. | v1.2 - corrected end-of-range code for IGOR Pro 3.01
  5. | v1.1 - Igor Pro version 2.05 fixes a bug with a or b at first or last point in wave,
  6. | and limits the range to a valid range if a or b are before or after the
  7. | first and last points in the wave.
  8. |
  9. | Given a tabulated function with monotonically increasing or decreasing
  10. | x values, AreaXY returns the trapezoidal area over an interval.
  11. |
  12. Function/D AreaXY(xwave,ywave,a,b)
  13.     Wave/D xwave,ywave                | xwave must be monotonic!
  14.     Variable/D a,b                        | limits of integration
  15.  
  16.     variable n= numpnts(xwave),tmp
  17.     variable/D increasing= xwave[n-1] > xwave[0],reversed=0
  18.     
  19.     if( increasing )
  20.         if( a > b )
  21.             tmp= a
  22.              a= b
  23.              b= tmp
  24.             reversed= 1
  25.         endif
  26.     else
  27.         if( a < b )
  28.             tmp= a
  29.              a= b
  30.              b= tmp
  31.             reversed= 1
  32.         endif
  33.     endif
  34.     variable/D pa=BinarySearch(xwave,a),pb=BinarySearch(xwave,b)
  35.     if( pa < 0 )
  36.         pa= 0
  37.         a= xwave[pa]
  38.     endif
  39.     if( pb < 0 )
  40.         pb= n-1
  41.         b= xwave[pb]
  42.     endif
  43.     pb += 1
  44.     | now we know that a and b are somewhere between xwave[pa] and xwave[pb] inclusive
  45.     variable/D f= (xwave[pa+1]-a)/(xwave[pa+1]-xwave[pa])
  46.     variable/D y= f*ywave[pa]+(1-f)*ywave[pa+1]                | interpolated y at a
  47.     variable/D area= (y+ywave[pa+1])*(xwave[pa+1]-a)/2        | area from a to next point
  48.     do
  49.         pa += 1
  50.         if( pa >= pb )
  51.             break;
  52.         endif
  53.         area += (ywave[pa]+ywave[pa+1])*(xwave[pa+1]-xwave[pa])/2
  54.     while(1)
  55.     if( pb != n )
  56.          f= (xwave[pb]-b)/(xwave[pb]-xwave[pb-1])
  57.         y= f*ywave[pb-1]+(1-f)*ywave[pb]                | interpolated y at b
  58.         area -= (y+ywave[pb])*(xwave[pb]-b)/2            | correct for area from b to last point
  59.     endif
  60.     if(reversed)
  61.         return -area
  62.     endif
  63.     return area
  64. End
  65.  
  66.